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ABSTRACT 

Spacetime Discontinuous Galerkin (DG) methods are used to solve hyperbolic PDEs describing wavelike physical 
phenomena. When the PDEs are nonlinear, the speed of propagation of the phenomena, called the wavespeed, at any 
point in the spacetime domain is computed as part of the solution. We give an advancing front algorithm to construct 
a simplicial mesh of the spacetime domain suitable for DG solutions. Given a simplicial mesh of a bounded linear 
or planar space domain M, we incrementally construct a mesh of the spacetime domain M x [0, oo) such that the 
solution can be computed in constant time per element. We add a patch of spacetime elements to the mesh at every 
step. The boundary of every patch is causal which means that the elements in the patch can be solved immediately 
and that the patches in the mesh are partially ordered by dependence. The elements in a single patch are coupled 
because they share implicit faces; however, the number of elements in each patch is bounded. The main contribution 
of this paper is sufficient constraints on the progress in time made by the algorithm at each step which guarantee that 
a new patch with causal boundary can be added to the mesh at every step even when the wavespeed is increasing 
discontinuously. Our algorithm adapts to the local gradation of the space mesh as well as the wavespeed that most 
constrains progress at each step. Previous algorithms have been restricted at each step by the maximum wavespeed 
throughout the entire spacetime domain. 

Keywords: mesh generation, unstructured meshes, advancing front, partial differential equations, 
discontinuous Galerkin, nonlinear hyperbolic PDE 



1. INTRODUCTION 

Simulation problems in mechanics consider the behav- 
ior of an object or region of space over time. Scien- 
tists and engineers use conservation laws and hyper- 
bolic partial differential equations (PDEs) to model 
transient, wavelike phenomena propagating over time 
through the domain of interest. Example applications 
are numerous, including, for instance, the equations of 
elastodynamics in seismic analysis and the Euler equa- 
tions for compressible gas dynamics. Closed-form so- 
lutions are typically unavailable for these problems, so 
analysts usually resort to numerical approximations. 

Finite element methods (FEM) are popular options 
for solving this class of problems. In the standard 
semi- discrete approach, a finite element mesh dis- 
cretizes space to generate a system of ordinary dif- 



ferential equations in time that is then solved by a 
time-marching integration scheme. Most semi-discrete 
methods impose a uniform time step size over the en- 
tire spatial domain, i.e., the time step does not adapt 
to the local gradation of the space mesh. Therefore, 
the resulting spacetime mesh consists of many more 
elements than required by physical causality. Hence, 
algorithms that use a nonuniform time step size can 
substantially improve computational efficiency. 

Spacetime discontinuous Galerkin (DG) methods have 
been proposed by Richter 8 , Lowrie et al. , and Yin 
et al. [11] for solving systems of nonlinear hyperbolic 
partial differential equations. Like traditional finite 
element methods, spacetime DG methods use basis 
polynomials to approximate the solution within each 
element; however, unlike traditional FEM methods, 
these basis polynomials have local support restricted 



to each element and the basis polynomials of adjacent 
elements do not have to agree on their common inter- 
section. This approach eliminates artificial coupling 
between adjacent elements when the mesh satisfies cer- 
tain causality constraints. (For further background on 
general discontinuous Galerkin methods, we refer the 
reader to Cockburn, Karniadakis, and Shu [3]-) 

Ungor and Sheffer flO] and Erickson et al. [4] devel- 
oped the first algorithm, called 'TentPitcher', to build 
graded spacetime meshes over arbitrary simplicially 
meshed spatial domains, suitable for spacetime DG so- 
lutions. Unlike most traditional approaches, the Tent- 
Pitcher algorithm does not impose a fixed global time 
step on the mesh, or even a local time step on small 
regions of the mesh. Rather, it produces a fully un- 
structured simplicial spacetime mesh, where the dura- 
tion of each spacetime element depends on the local 
feature size and quality of the underlying space mesh. 

Efficient spacetime meshing relies on the notion of the 
domain of influence and the domain of dependence of 
an event. Imagine dropping a pebble into a pond — 
circular waves propagate outwards from the point of 
impact. The frontier of expanding waves sweeps out 
a cone in spacetime called the domain of influence of 
the event. The radius of the domain of influence at 
any time is the radius of the circular disc consisting 
of all points on the surface where the initial wave has 
arrived. The domains of influence and dependence can 
be approximated by right circular cones with common 
apex P (Figure [T]). The symmetric double cone rep- 
resenting the domains of influence and dependence at 
points P in spacetime can be described by a scalar fleld 
uj where uj{P) — dr/dt, the wavespeed at P, specifies 
how quickly the radius r of domains of influence and 
dependence of P grows as a function of time. Smaller 
values of aj(P), i.e., steeper cones, correspond to slower 
wavespeeds. The wavespeed u}{P) at a point in space- 
time is, in general, part of the solution of the PDE 
at that point. The slope of the cones of influence and 
dependence of P, denoted by cr(P), is the reciprocal of 
the wavespeed — larger slopes mean steeper cones and 
therefore slower wavespeeds, and smaller slopes mean 
shallower cones and faster wavespeeds. 

Given a simplicial mesh of some bounded domain 
M C R'*, the Tent Pitcher algorithm incrementally 
constructs a simplicial mesh of the spacetime domain 
using an advancing front method. The spacetime do- 
main is the subset M x [0,00) C R'^'''^, a subset of 
Euclidean space one dimension higher. The algorithm 
progresses by adding simplices to the evolving mesh 
in small patches by moving a vertex of the front for- 
ward in time. The inflow and outflow boundaries of 
each patch (Figure [2]) are causal by construction, i.e., 
each boundary facet F separates the cone of influence 
from the cone of dependence of any point on F (Fig- 



ure [TJ. Equivalently, for every point P on F we have 
||VP|| < l/aj(P) = cr(P). If the outflow boundaries of 
a patch are causal, every point in the patch depends 
only on other points in the patch or points of inflow el- 
ements adjacent to the inflow boundaries of the patch. 
Therefore, the solution within the patch can be com- 
puted as soon as the patch is created, given only the 
inflow data from adjacent inflow elements. The ele- 
ments within a patch are causally dependent on each 
other and must be solved as a coupled system. Pro- 
vided the space mesh has constant degree, each patch 
contains only a constant number of elements and can 
therefore be solved in constant time. Therefore, the 
computation time required to compute the numerical 
solution is linear in the number of spacetime elements. 
Patches with no causal relationship can be solved in- 
dependently. To minimize undesirable numerical dis- 
sipation and the number of patches, we would like the 
boundary facets of each patch to be as close as possible 
to the causality constraint without violating it. 

The causality constraint limits the progress in time 
at each step, i.e., the height of each tentpole is con- 
strained. For spatial domains of dimension d > 2, it is 
not trivial to guarantee that the advancing front algo- 
rithm can always make progress. We require that for 
any target time value T the algorithm will compute 
a mesh of the spacetime volume M x [0, T] and the 
solution everywhere in this volume in finitely many 
steps. The target time T is not known a priori be- 
cause it depends on the evolving physics. The original 
Tent Pitcher algorithm proposed by Ungor and Shef- 
fer [TD] applied to one- and two-dimensional space do- 
mains. The algorithm could guarantee progress only 
if the input triangulation contained only angles less 
than 90 degrees and if the wavespeed did not increase 
or increased smoothly. Erickson et al. [4] extended 
Tent Pitcher to arbitrary spatial domains in any di- 
mensions by imposing additional constraints, called 
progress constraints. The progress constraint applied 
to a single simplex on the front limits the amount of 
progress in time when some vertex of the simplex is 
pitched. The progress constraint is a function of the 
shape of the simplex. The geometric constraints that 
limit the height of each tentpole are called cone con- 
straints. 

All the results so far have applied to the case where 
the wavespeed at a given point is either constant, de- 
creasing, or increasing smoothly as a Lipschitz func- 
tion. (See Alper Ungor's PhD thesis [9] for the details.) 
When the wavespeed changes, the previous algorithms 
take the fastest that the wavespeed can ever be and 
use that as a conservative upper bound on the wave- 
speed at any time. One would like an algorithm that 
adapts to increasing wavespeeds so that fewer space- 
time elements, and therefore less computation time, 
are required to mesh a given volume. 



In this paper, we give an advancing front algorithm 
to construct a spacetime mesh over an arbitrary hn- 
ear or planar space mesh (d < 2). Our algorithm 
extends TentPitcher to the case when the wavespeed 
can be an arbitrary scalar field over the spacetime do- 
main. In particular, our algorithm guarantees finite 
positive progress at each step even when the wave- 
speed at a given point increases discontinuously and 
unpredictably over time. 

The main contributions of this paper are twofold. We 
give a novel characterization of fronts that are al- 
ways guaranteed to progress, which we call progressive 
fronts, and give a lower bound on the progress guaran- 
tee at each step which depends only on the local size of 
the mesh and the wavespeed that most constrains the 
duration of the current patch. The minimum progress 
guarantee at any step is a positive quantity bounded 
away from zero, so the front is guaranteed to progress 
past any target time in a finite number of steps. The 
second contribution of this paper is to give geomet- 
ric constraints on the front at any step that guarantee 
that the front can progress in the next step and so on 
inductively at every step. The geometric constraints 
are simple to express and to compute. Intuitively, the 
geometric constraints that apply at any given itera- 
tion of the algorithm are predicted by looking ahead 
at the next iteration of the algorithm. We also give an 
efficient algorithm to maximize the progress at every 
step subject to these constraints. The novelty of our 
characterization of progressive fronts and of our algo- 
rithm is that we resolve the following conundrum. The 
progress of the front at each step i is limited by the 
progress constraint that must be satisfied by the next 
front at step i -I- 1. However, we do not know what is 
the next front unless we know how much progress is 
possible at step i. 

The paper by Erickson et al. 2] contains an error in 
the statement of the causality constraint when obtuse 
triangles are involved; therefore, their proof of correct- 
ness is incomplete because it omits the obtuse angle 
case. While their proof can be fixed, we prefer our 
new algorithm, which is provably correct even when 
the wavespeed is constant or does not increase. Our 
new progress constraints are potentially weaker than 
those of Erickson et al. [4] . 

Our algorithm is the first algorithm to build space- 
time meshes over arbitrary planar triangulated spa- 
tial domains suitable for solving nonlinear hyperbolic 
PDEs, where the wavespeed at any point in spacetime 
depends on the solution and cannot be computed in 
advance. Moreover, the solution can change discontin- 
uously, for instance when a shock propagates through 
the domain. 

The input to our advancing front algorithm is a simpli- 
cially meshed bounded domain M C R** where d < 2 




Figure 1: A causal face separates the cones of influence 
and dependence at every point on the face. 
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Figure 2: A vertical cross-section of a patch of tetrahe- 
dra; the inflow and outflow faces are causal. 

and the initial conditions of a nonlinear hyperbolic 
PDE. The space mesh describes the situation at time 
equal to zero, specifically, the slope at every point in M 
at time zero. We allow more general initial conditions 
but we will postpone a description of those conditions 
until later sections. Our meshing algorithm is an ad- 
vancing front procedure which alternately constructs 
a new patch of elements and invokes a spacetime DG 
finite element method to compute the solution within 
that patch. At every iteration, the front is the graph of 
a continuous piecewise linear time function t : M ^ M.. 
The front t is linear within every simplex of M and 
llVt(p)|| < o"(p) for every point p £ M. The front is a 
terrain whose facets correspond to simplices in the un- 
derlying space mesh. Each facet of the front coincides 
with the outflow face of a patch in the past and the in- 
flow face of a patch in the future. We say that a front 
is causal if every simplex of the front is causal. To ad- 
vance the front t, the algorithm chooses an arbitrary 
vertex P — {p, t(p)) from the front and lifts it to a new 
point P' = {p,t'{p)) where t'{p) > t{p) and for every 
other vertex q we have t'{q) — t{q). The spacetime 
volume between the new front t' and the old front t is 
called a tent. The tent is meshed with simplices shar- 
ing the edge {P, P') called the tentpole. The height of 
the tentpole is the duration t'{p) — t{p). Consider a 
planar space mesh M. For each triangle pqr incident 
on p, the tetrahedron P' PQR belongs to the patch. 
The outflow face P'QR and the inflow face PQR are 
causal boundaries. The triangles P'PQ and P' PR are 
implicit faces. Since the implicit faces are vertical they 



are not causal boundaries and so elements within the 
patch are coupled. The elements below the front t 
whose outflow faces intersect any of the inflow faces of 
the new patch are inflow elements. We pass the newly 
constructed patch along with all its inflow elements to 
a DG solver. The DG solver returns as part of the 
solution the slope at every point on every outflow face 
of the patch. The new front t' and the output of the 
DG solver are the input to the next iteration of the 
algorithm. 

Since we are interested in causal fronts only, hence- 
forth it is implicit that every front considered is causal. 

We assume that the slope at any point P is bounded 
by the minimum and the maximum slopes anywhere 
in the cone of dependence of P. Hence, given a front t 
and a point P in the future, the slope at P is no smaller 
than the slope at Q for every point Q on the front t 
such that P is in the cone of influence of Q. 

It can be computationally very expensive to determine 
the shallowest cone of influence that contains a given 
point P. In particular, the shallowest cone of influence 
containing P may correspond to a nonlocal point Q, 
one arbitrarily distant from P. To compute this non- 
local cone constraint eSiciently, we use a standard hi- 
erarchical decomposition, called a bounding cone hi- 
erarchy, of the space domain. The elements in the 
hierarchy correspond to subsets of the space domain. 
For each element of the hierarchy, we compute the 
minimum slope within the corresponding subset of the 
space domain. The smallest element in the hierarchy 
is a single simplex. In order to determine the strictest 
cone constraint that applies locally, we traverse the 
hierarchy until we determine the simplex with min- 
imum slope whose cone of influence contains P. In 
practice, we expect that our algorithm has to exam- 
ine only a small subset of the hierarchy. In the worst 
case, the algorithm has to examine every simplex of the 
front but in that case the algorithm will be at most a 
constant factor slower than one that does not use a 
bounding cone hierarchy. When a patch is solved, the 
bounding cones are updated with the new slopes by 
traversing a path from a leaf to the root of the hier- 
archy. This hierarchical approximation technique has 
been applied very successfully to numerous simulation 
problems, such as the Barnes-Hut divide-and-conquer 
method [2] for A'^-body simulations, as well as to colli- 
sion detection in computer graphics and robot motion 
planning ^ and for indexing multi-dimensional data 
in geographic information systems [Sj. 

1.1 Notation 

We use lowercase letters like p, q, r to denote points 
in space and uppercase letters like P, Q, R to denote 
points in spacetime. A front £ is a piecewise linear 



function i : M — > R. For a simplex (of any dimensions) 
T of M, let t\^ denote the time function t restricted to 
T and extended to the affine hull of r; in other words, 
t\^ is a linear function that coincides with t for every 
point of T. Let ii : M — > R denote the front after 
the ith step of the algorithm; to is the initial front. 
For every i, the front ti is a terrain whose facets are 
the simplices of M. In other words, ti is a piecewise 
linear function such that for every simplex r of M, the 
functions ti and ti\^ coincide at the vertices of r. 

For a time function £ : Af ^ R we denote the gradient 
of f by V £. A local minimum of the front £ is a vertex p 
such that t{p) < t{q) for every vertex q that is a neigh- 
bor of p. When the current front t is clear from the 
context, for every point p £ M we use P to denote the 
corresponding point on the front, i.e., P — {p,t{p)). 

For a point P in spacetime, we use o{P) to de- 
note the reciprocal of the wavespeed at P. Let 
(Tmin denote minpgMx[o,c3o){o-(P)} and amax denote 
maxpgMx[o,oo){o"(f')}- We assume that < (Tmin < 
o"max < oo. For a simplex r in spacetime, we use (j(r) 
to denote the minimum of fj{P) over all points P in r. 

We say that a front t' is obtained by advancing a ver- 
tex p of M by 5t > if t'{p) = t{p) -\- 5t and for every 
other vertex q ^ p we have t'{q) = t{q). For any front 
t, vertex p, and real St > 0, let t' = next{t,p, St) de- 
note the front obtained from t by advancing p by St. 

1.2 Problem statement 

The input to our problem is the initial front to and the 
initial conditions of the PDE. We want an advancing 
front algorithm such that for every T G R-'^ there 
exists a flnite integer fc > such that the front tk after 
the fcth iteration of the algorithm satisfies tk > T. 

We say that a front t is valid if there exists a posi- 
tive real S bounded away from zero such that for ev- 
ery T G R-" there exists a sequence of fronts t, ti, 
t2, . . ., ife where tk > T, each front in the sequence 
obtained from the previous front by advancing some 
vertex by S. What makes the definition of a valid front 
nontrivial is the requirement that all fronts be causal. 
The main difficulty in characterizing valid fronts arises 
when the wavespeed at a given point in the space do- 
main increases discontinuously and unpredictably over 
time. 

Our solution We define progressive fronts and prove 
that if a front is progressive then it is valid. We give 
an algorithm that given any progressive front ti con- 
structs a next front ti+i such that ii+i is progressive. 
The volume between ti and ti+i is partitioned into sim- 
plices. The next front ti+i is obtained by lifting a local 
minimum of ti by a positive amount bounded away 
from zero. The algorithm can easily be parallelized 



to solve several patches asynchronously by lifting any 
independent set of vertices in parallel. Whenever the 
algorithm chooses to lift a local minimum, it is guar- 
anteed to be able to lift it by at least Tmin > which is 
a function of the input and bounded away from zero. 

2. ONE-DIMENSIONAL SPACE 
DOMAINS 

We begin by describing our algorithm to construct 
spacetime meshes over one-dimensional space do- 
mains. Even this simple case captures all but one as- 
pect of the complexity of guaranteeing causality when 
wavespeeds are changing. 

The space domain M is a closed interval of the real 
line. The input space mesh is a subdivision of this 
interval into segments. Let V{M) denote the set of 
vertices of the space mesh M . The initial front to cor- 
responds to io(p) = for every vertex p of the space 
mesh, but more generally, any (causal) front can be 
the initial front. Let Wmin denote the minimum length 
of any segment in the space mesh. Let (Tmin denote the 
minimum slope (t{P) over every point P in the space- 
time domain M x [0,oo). Let Tmin denote (T,ninWmin. 

In iteration i + 1 of our advancing front algorithm 
(i > 0), we advance a single vertex p, where p is a 
local minimum of the current front ti, to get the new 
front fi+i, i.e., ti+i = next(t,p, Jt). More generally, we 
can advance any vertex or an independent set of ver- 
tices, not necessarily local minima, forward in time. 
The value of ti+i{p) is bounded from above by the 
requirement that ti+\ be causal. 

Let AB be an arbitrary segment of the front fi+i. 
Without loss of generality, assume ti+i{a) < ti+i{b). 
Then, AB is causal if and only if the gradient of the 
time function ti+i restricted to ab is at most the slope 
a{AB), i.e., if and only if 



Theorem 1. Let ti be a front and letp be an arbitrary 
local minimum ofti. Then, for every St £ [0, Tmin] the 
front ti-f-i = next(ti,p, St) is causal. 

Proof. Only the segments of the front incident on P 
advance along with p. Consider an arbitrary segment 
pq incident on p. Let t and t' denote til^^ and ti+i\^^ 
respectively. We have t{p) + St < t{p) + (Tminifmin < 
t{q) + \pq\a{P'Q) because p is a local minimum, uimin < 
\pq\, and amin < a{P'Q). Therefore, the segment P'Q 
is causal. Since this is true of an arbitrary segment 
on the front t' , we have proved that that the front 
ti+i = next(ti,p, 5t) is causal. □ 



Theorem 2. For any i > 0, if the front ti is causal 
then ti is valid. 

Proof. Consider step z + 1 of the algorithm. By The- 
orem [1] the front ti+i such that ti+i{p) £ [0,Tmin] 
is causal. Therefore, we have shown that if ti is 
causal then there is a front ti+i — next(ti,p, Tmin) 
such that ti+i is causal. Note that X]pGV{M) = 
rmin + X]pgv{M) (p) • By Inductlon on i, and because 
o"max is finite and A'l is bounded, there exists a finite 
k > i such that the front tk satisfies 

J2 ik{p) > diam(Af)(JmaxT 

pev(M) 

for any real T. Since tk is causal 

( max tfc(p) ) < diam(M)amax ( min tfc(p) ) . 
Therefore, minpgv'(M) ife(p) > T and so ti is valid. □ 

2.1 Being greedy at every step 

We would like to maximize the progress at each step 
in a greedy fashion, i.e., given a front ti we would 
like to maximize ti+i{p) where ti+i = next{ti,p, St) 
subject to the constraint that ti+i is causal. By The- 
orem[2] we can have ti+i(p) > ti{p) + Tmin- However, 
it may be possible to make further progress by setting 
ti+i{p) higher, especially if each segment PQ incident 
on p each satisfies progress constraint [cTprcv] for some 
fprcv < ct{P'Q) at the end of the previous iteration. 

For a fixed segment pq incident on p let T^^^ denote 
sup{r : P'Q is causal where P' = (p,T)}. To max- 
imizing the progress at step i + 1, we would like to 
compute TsJp'. The segment P'Q is causal if and only 
if the slope of P'Q is less than or equal to the slope of 
the cone of influence from every point on the front that 
intersects P'Q. A cone of influence intersects P'Q if 
and only if the cone intersects the tentpole PP' . In 
general, a cone of influence from arbitrarily far away 
can intersect the tentpole at p. See Figure [S] This 
is not the case when the wavespeed everywhere is the 
same. Therefore, in general, Tl^^ could be determined 
by a cone of influence of a point arbitrarily distant 
from p. 

Partition the front into two subsets of points: 

(i) points in the star of P ("local" points), and 

(ii) points everywhere else on the front ("remote" 
points). Corresponding to each subset we have 
two disjoint subsets of cones of influence — Ciocai and 
Crcmote respectively. Each subset of cones limits the 
new time value of p and so the final time value is the 
smaller of the two values for each of Ciocai and Crcmote 
taken separately. 




Figure 3: Top to bottom: a sequence of tent pitch- 
ing steps in IDxTime. IVIaximizing the height of each 
tentpole while staying below every cone of influence can 
require examining remote cones arbitrarily far away. 



Consider the subset Ciocai- Let aiocai denote the small- 
est slope among all cones of influence in Ciocai The 
segment P'Q is causal only if its slope is less than or 
equal to criocai- Let Tiocai be the maximum time value 
of P' for which the slope of P'Q is less than or equal 
to (Jiocai- The maximum Tiocai exists because the set 
of feasible values is closed and therefore compact. To 
compute Tiocai we substitute criocai in the condition for 
causality of P'Q (Equation [l| . 

Next consider the subset Cremote- The front ti is 
strictly below every cone in demote because ti is causal. 
The segment P'Q is causal only if it is also strictly be- 
low every cone in Cremote- Given a cone C € Cremote, C 
intersects P'Q if and only if C intersects the tentpole 
PP' . Let Tremote deuote the smallest time value T for 
which the tentpole PP' where P' — (p,T) intersects 
exactly one cone in Cremote- The segment P'Q is causal 
only if T < Tremote- Note that the upper bound on T 
imposed by remote cones is a strict inequality. 

Therefore, the progress ti+i{p) — ti{p) at step i + 1 
is limited because Tg+,^ = max{Tiocai, Tremote}- To 
maximize the progress at the current step, we choose 
ti+i{p) equal to T^'r^p^ minus the machine precision 77, 
or ti{p) + Tmin, whichever is larger- 
Computing Tremote eXHCtly Computing Tremote 

is equivalent to answering a ray shooting query in the 
arrangement of the cones in Cremote- We use a bound- 
ing cone hierarchy TC obtained from a hierarchical de- 
composition of the space domain to efficiently answer 
the ray shooting query- The hierarchical decompo- 
sition of the space domain induces a corresponding 
hierarchical decomposition of every front. For each el- 
ement of this hierarchy, we store a right circular cone 
that bounds the cone of influence of every point of the 
corresponding subset of the front. To answer the ray 
shooting query, we traverse the cone hierarchy from 
top to bottom starting at the root. At every stage, 
we store a subset C of bounding cones such that every 
cone in Cremote is contained in some cone in the subset 
C. The cones in C are stored in a priority queue in 
non-decreasing order of the time value at which the 
vertical ray at P intersects each cone. Initially, C con- 
sists solely of the cone at the root of the hierarchy. 
At every stage, if the cone in C that has the earliest 
intersection does not come from a leaf in the hierar- 
chy then we replace it in the priority queue with its 
children. Continuing in this fashion, we eventually de- 
termine the single facet of the front such that the cone 
of influence from some point on this facet is intersected 
first by the vertical ray at P. The time coordinate of 
the point of intersection is Tremote, the answer to the 
ray shooting query. 

If the hierarchy is balanced its depth is 0(log m) where 
m is the number of simplices in the space mesh. In 



ID X Time, we observed empirically that on average 
only a few nodes in the cone hierarchy were examined 
by this algorithm to determine the most constraining 
cone of influence. 

Approximating Tremote Since we know a range of 
values [ti(p)+Tmiii, Tiocai] that contains Trcmotc, we can 
approximate Trcmotc up to any desired numerical accu- 
racy by performing a binary search in this interval. At 
every iteration, we speculatively lift P to the midpoint 
of the current search interval. Let P" be the specula- 
tive top of the tentpole at P. We query the cones of 
influence in Cremote to determine the minimum slope 
fremote among all cones that intersect PP" . If the 
maximum slope of the outflow faces incident on P" is 
less than (Jrcmotc then we can continue searching in the 
top half of the current interval; otherwise, the binary 
search continues in the bottom half of the current in- 
terval. The search terminates when the search interval 
is smaller than our desired accuracy. A bounding cone 
hierarchy helps in the same manner as before to de- 
termine the minimum slope among all cones in demote 
that intersect PP". 

Theorem 3. Given a simplicial mesh M of a bounded 
real interval where Wmin is the minimum length of 
a simplex of M and Umin is the minimum slope 
anywhere tn M x [0,cx3) our algorithm constructs 
a simplicial mesh of M x [0, T] consisting of at 
most j" ^ '^"'"'^"'l^^ """""^ spacetime elements for every 
real T > 0. 

Proof. In Theorem[T] we have shown that the height of 
each tentpole constructed by the algorithm is at least 
Tmin = fmin^^min. By Thcorcm [J] after constructing 
at most k < j''^ patches, the entire front 

tk is past the target time T. Since each patch consists 
of at most two elements, the theorem follows. □ 

We have shown that every causal front in ID x Time is 
valid. In higher dimensions, additional progress con- 
straints are necessary. 

3. PLANAR SPACE DOMAINS 

In this section, we describe our algorithm for d — 2, 
i.e., for a triangulated planar space domain M C K^. 

For planar domains, we encounter nontrivial progress 
constraints that are necessary to guarantee sufBcient 
progress at each step, i.e., to guarantee that the height 
of the tentpole constructed at every step is positive 
and bounded away from zero. In the absence of such 
constraints, it was shown by Ungor and Sheffer [lOj . 
and by Erickson et al. [4] that if the space mesh con- 
tains an obtuse or a right triangle then Tent Pitcher 



will eventually construct a front such that no further 
progress is possible while maintaining causality. Er- 
ickson et al. [3] derived additional progress constraints 
that were sufficient to guarantee progress, even in the 
presence of obtuse angles, however only by assuming 
the minimum slope occurs everywhere in spacetime. 
In this section, we show how to relax these progress 
constraints so that they adapt to the slope of the 
most constraining cone of influence at every step. Our 
progress constraint is a function of the slope encoun- 
tered locally in the next step of the algorithm, which 
may be substantially less constraining than the glob- 
ally minimum slope. 

Fix a real parameter e G (O, |] . The space domain M 
is a triangulation of a bounded subset of the plane . 
Let Wmin denote the minimum width of any triangle of 
the space mesh. Let Omin denote the minimum ij{P) 
over every point P in the spacetime domain A/ x [0, oo) . 

Let Tmin denote eCTminWmin. 

Definition 1 (Progress constraint a). Let PQR be an 

arbitrary triangle of a front t. Without loss of gener- 
ality, assume t(p) < t{q) < t(r). We say that the tri- 
angle PQR satisfies progress constraint a if and only 

where (j)p = max {sin Zprg, sin Zpgr}. Note that < 

4>p < 1- 

Suppose the lowest vertex p is being advanced. As 
long as p is the lowest vertex of Apqr, the progress 
constraint limits ||V t\^^\\ but ||V t\^^\\ is unchanged 
by lifting p. When t{p) > t{q), the new lowest vertex 
is q, so the progress constraint limits ||V t\^^\\. (We 
can interpret the progress constraint inductively as a 
causality constraint on the 1-dimensional facet pr op- 
posite q where the relevant slope is (1 — e)a<j)q.) 

Definition 2 (Progressive). Let t be a front and let 
pqr be a given triangle. Without loss of generality, 
assume t{p) < t{q) < t{r). We say that the triangle 
PQR is progressive if and only if both of the following 
conditions are satisfied by P'QR where P' = (p, t{p) -\- 
St) for every St G [0, Tmin]: 

1. P'QR is causal, and 

2. P'QR satisfies progress constraint a{P'Q'R) 
where Q' = {q,t(q) -\- T,m„). 

We say that a front t is progressive if every triangle on 
the front is progressive. Note that every progressive 
triangle or front is also causal. 



3.1 A new advancing front algorithm 

We are now ready to describe iteration i + 1 of our 
advancing front algorithm for i > 0. Advance a sin- 
gle vertex p by a positive amount, where p is any lo- 
cal minimum of the current front ti, to get the new 
front ti+i such that for every triangle pqr incident on 
p the corresponding triangle on the new front ti+i is 
progressive. In the parallel setting, advance any in- 
dependent set of local minima forward in time, each 
subject to the above constraint. The value of ti+i{p) 
is constrained from above separately for each of the 
simplices incident on p. The final value chosen by the 
algorithm must satisfy the constraints for each such 
triangle. Therefore, it is sufficient to consider each 
triangle pqr incident on p separately while deriving 
the causality and progress constraints that apply while 
pitching p. 

Next, we derive simple formulae for the causality and 
progress constraints for a given triangle pqr when p is 
being pitched. Let t and t' denote ti\^^^ and U+il^^^ 
respectively. 

Let ftqr denote the unit vector normal to qr such 
that nqr-(p—q) > 0. Let Vqr be the unit vector parallel 
to qr such that Vqr ■ {f— q) > 0. Then, {ugr, Vqr} form 
a basis for the vector space R^. Let n^p denote the 
unit vector normal to pr such that firp ■ [q — p) > 0. 
Let Vrp be the unit vector parallel to rp such that 
fp ■ (p — r) > 0. Then, {nrp,Vrp} form another basis 
for the vector space R^. 

The gradient vector V t' can be written as 



where 



Vf' = (Vt' ■ nqr)nqr "f V t' j 



Lifting p does not change the gradient of the time func- 
tion restricted to the opposite edge, so Vt' • Vqr = 
Vt ■ Vqr, i.e., V t'\^^ = V t\^^.. Since q is the lowest 
vertex of qr, we have V t' ■ Vqr = V t • Vqr > 0. 



Also, 



where 



V t' — {V t' ■ nrp)nrp 4- V t'j 



The vectors n^r and n,.p are related by a rotation 
around the origin by angle 6. Since < 6^ < tt we 
have cos 6' = n,,- ■ rtrp and sin^ = \/l — {uqr ■ Urp)^. 
Hence, 



l^*Xpll = \\^t\qr\\ COS e + {V t' -nqr) Sine 

= II V t\^J{nqr ■ Urp) 



+ {Vt' ■ nqr)^l ~ {f. 



(2) 






(a) (b) (c) 

Figure 4: Triangle pqr where t{p) < t{q) < t{r) 

Deriving the causality constraint Let u be the 

orthogonal projection of p onto line qr. Since lifting p 
does not change the time function restricted to qr, we 
have t'lg^ = t\^^. The scalar product \7 1' ■ fiqr can be 
written as 

t\p)-t{u) 



V t' ■ fiqr = 



\up\ 



Since q is the lowest vertex of qr and since PQR is 
progressive, we have < V t' • Vqr ~ V t ■ Vqr < 
{I - £)a{P'QR) < a{P'QR). Therefore, ||Vt'|| < 
a{P'QR) if and only if 



t'{p)-t{u) 



\up\ 



< Ja{P'QRY 



\V t\ IP 



(3) 



Deriving the progress constraint Let (Tprog de- 
note a{P'Q'R) where P' = {p,t'{p)) and Q' = 
(g, -j-Tmin). By Equationg] the triangle P'QR sat- 
isfies the progress constraint ||V < (1 — e)(7prog</'q 
if and only if 

^ ^ (1 - e)c^prog0q - ||V t| • ^rp) 

V t ■ Uqr < ^ 

^/l — {fiqr ■ firpY 

Therefore, the progress constraint is 



t'{p)-t{u) 

I 



■{l^£)0{P'Q'R)(j}q 



^i-{a,r-ar„)^ " 



(4) 



3.2 Proof of correctness 



In this section, we prove the correctness of our algo- 
rithm, i.e., that every front constructed by the algo- 
rithm is valid. 

Theorem 4. // a front ti is progressive, then for any 
local minimum vertex p and for every St G [0, T-min\ the 
front fi+i = next(ti,p,St) is causal. 



Proof. Since only the triangles of the front incident 
on P advance along with p, we can restrict our atten- 
tion to an arbitrary triangle pqr incident on p. Let t 
and t' denote ti\^^^ and ti+il^^^ respectively. Let u 
be the orthogonal projection of p onto line qr. 



Consider the causality constraint (Equation |3]). We 
will consider two cases separately: (i) t(u) > t{q), and 
(ii) t{u) < t{q). 

Case 1: t(u) > t(q) > t(p) See Figure |li;b)-(c). In 

this case, we have 

t'{p)^t{p) + St 

< t{u) + St 

< t{u) + ea{P'QR)\up\ 

because \up\ > Wmin and a{P'QR) > Umin- Since < 
e < i we have e < ^1 — (1 — e)^. Therefore, 

t'{p) < t{u) + \up\^l-{l-e)2a{P'QR) 

= t{u) + \up\^a{P'QRy' - (1 - e)2cr2(P'Q7?) 



which is precisely the causality constraint of Equa- 
tion O The last inequality follows because PQR is 
progressive, hence |1V < [1 — e)a{P'QR)(j>p < 

{l-e)a{P'QR). 

Case 2: t(u) < t(q) See Figure i^a). Let /3 = 
\uq\/\up\. Since \uq\ 7^ 0, we have 

t'{p) - t{u) t'{p) - t{q) , t{q) - t{u) \uq\ 



\up\ 



\up\ \uq\ \up\ 

t'{p)-t{q) 



\up\ 



(5) 



Using Equation [S] the causality constraint (Equa- 
tion |3]) can be rewritten as 



t'{p)-t{q) 



\up\ 



<Ja{P'QRr^\\Vt\\\^ 



(6) 



Since ti is progressive, we have ||V t\^^\\ < (1 — 
e)<j{P'QR)(j>p. Substituting this upper bound on 
||V t\^^\\ into Equation[Sl we obtain the following con- 
straint: 

^^^^ < '^iP'QR) (Vi-(i-)Vi) 

-a{P'QR){l-e)P<l>p (7) 

which implies the causality constraint of Equation |6] 
Now, 



t'{p)-t{q) ^ t'{p)-t{p) 



\up\ 



< 



\up\ 

£^ i n ^mi n 



\up\ 

< eo-min. 



Since (Tmin < a{P'QR), Equation [7] is satisfied if 
or equivalently 



We have 



(e + (1 - e)[3(l,pf + {l-ef4>l = l + 2e(l-e) I 



-1) 



We have (j>p = sin Zpgr = \up\/\pq\ > \up\/\pr\ — 
sin Zprq and /3 — |ug|/|up|. Since \uq\ < \pq\, we have 
Ptjjqr < 1- Therefore, Equation[7]is satisfied. □ 

Theorem 5. // a front t is progressive, then for any 
local minimum vertex p and for every St G [0, Tmin] the 
front t' — next{t,p,St) is progressive. 

Proof. Since only the triangles of the front incident 
on P advance along with p, we can restrict our atten- 
tion to an arbitrary triangle pqr incident on p. Let t 
and t' denote ti\^^^ and ti+i\^^^. respectively. Let u 
be the orthogonal projection of p onto line qr. Let 
CTprog denote a(P'Q'R) where P' — {p, t{p) + St) and 
Q' = {q,t{q) + Tn,in). 

We separate the analysis into three cases depending 
on which, if any, of the angles Zpqr and Zprq of Apqr 
is obtuse. 

Case 1: Both Zpqr and Zprq are non-obtuse. 
See Figure IHb). In this case, we have t{u) > 
til) ^ t{p) and fiqr ■ firp < 0. Let Q = 
y^l — {fiqr ■ fhrpY- Hcnce, siu Zqrp = a = \up\/\pr\ 
and Uqr ■ n,.p = — Vl — — —\ur\/\pr\. Also, cj)q = 
max {sin Zqrp, sin Zqpr}. 

Therefore, the progress constraint of Equation |4] can 
be rewritten as follows: 

t'(p)-t(u) ^ (1 c)a(P'C' R) '"=^''{°'"^'?''P'°'"^'?P''} 

— V / V ^ / sin Zoru 



\up\ 



We have 



+ f^l|Vtl II 



(8) 



t'{p) - t{u) t'{p) - t{p) gg-minWmin 

^ II ^ II ^ g^min. 



\up\ 



\up\ 



\up\ 



Since e < i, we have e < 1 — e; also, Umin < a{P'Q' R); 



hence, 



£cr„,i„ < (1 - e)a{P'Q'R) 



< (1 - s)a{P Q R) 

< (1 - ,)^(p'Q'^) max{sinZgrp, sin Zgpr} 



sin Zqrp 



Therefore, the progress constraint of Equation[8]is sat- 
isfied. 

Case 2: Zpqr is obtuse. See Figure |4ja). In this 
case, we have t{u) < t{q) and fiqr ■ Urp < 0. Let a = 
■y/l — {Hqr ■ UrpY- Hence, sin Zqrp — a — |Mp|/|pr| 
and riqr ■ firp — — Vl — ~ —\ur\/\pr\. 

Let /3 = \uq\/\up\. Since \uq\ 7^ 0, we have 

t'jp) - t{u) _ t'{p)^t{q) ^ t{q)^t{u) \uq\ 
\up\ \up\ \uq\ \up\ 

_ t'{p) - t{q) 



\up\ 



+m^t\\\ 



Therefore, the progress constraint of Equation |4] can 
be rewritten as foUows: 



As before, we have 



t'{p) - t{q) t'(p) - t(p) ea-minWmin 



\up\ 



\up\ 



\up\ 



Since e < |, we have e < 1 — e; also, Umin < a{P'Q' R); 



hence, 



£cr„,i„ < (1 - e)a{P'Q'R) 



< (1 - e)a{P Q R) 

< (1 - ,),(^p'Q'R) ^'^-{^^^^irP,^^^^1Pr} 



sin Zqrp 



+ P 



\up\ 



llVti 



t'{p)-t{q) 



We have 



< (l-e)a(PVi?)" 



ax{sin Zgrp,sin Zqpr} 
sin /Lqrp 



(9) 



t'{p)~t{q) t'{p)-t{p) eO-minl«min , 

< : : < : : < eOmin 



\up\ 



\up\ 



\up\ 



Since e < ^, we have e < 1 — e; also, Umin < a{P'Q' R); 
hence, 

ecr^in < (1 - e)a(P'Q'R) 

< (1 - e)(j{P Q R) : — 

sm Zgrp 

< (1 - c)aiP'Q'R) "^^^^^^"^ ^'^'"P' ^^"^ ^QPr} 
~ sin Zqrp 



Therefore, the progress constraint of Equation O is 
satisfied. The last inequality follows because j3 = 

\uq\/\up\ < |iir|/|Mp|. 

Case 3: Zprq is obtuse. See Figure |3Ic). In this 
case, we have t{u) > t{r) > t{q) > t(p) and riqr ■ rirp > 
0. Let a — \/l — {fiqr ■ rirp)^. Hence, sin Zqrp — a — 
\up\/\pr\ and fiqr ■ firp = \/l — — \ur\/\pr\. 

Let /3 — \uq\/\up\. Since |itg| 7^ 0, we have 

t'{p) - t{u) _ t'{p) - t{q) , t(q) - t{u) \uq\ 



\up\ 



\up\ ^ \uq\ \up\ 
t'{p)-t{q) 



\up\ 



mtirW 



Therefore, the progress constraint of Equation |4] can 
be rewritten as follows: 



t'{v)-t{i) 



f\f r»\ max{sin Zgrp.siii zlqpr} 



I < {l^e)a{P'Q'R)i . , 



(10) 



Therefore, the progress constraint of Equation [10] is 
satisfied. The last inequality follows because P = 

\uq\/\up\ > \ur\/\up\. □ 

Theorem 6. For any i > 0, if the front tt is progres- 
sive then ti is valid. 

The proof is almost identical to that of Theorem (2] 



3.3 Being greedy 

We would like to maximize the progress at each step in 
a greedy fashion, i.e., given a front ti we would like to 
maximize ti+i{p) where ti+i — next{ti,p, 5t) subject 
to the constraint that ti+i is causal. For a fixed trian- 
gle ppiP2- ■ -Pd incident on p let T^^p denote sup {T : 
P'QR is causal and progressive, where P' = (p, T) and 
Pi = {pi,ti{pi) + Tmin}. To maximizing the progress 
at step i + 1, we would like to compute T^i^p . Sim- 
ilar to the ID X Time case, partition the set of cones 
of influence from points on the front ti into local and 
remote subsets. Let criocai denote the smallest slope 
among all local cones of influence. The triangle P'QR 
is causal only if its slope is less than or equal to criocai. 
Let Tiocai be the maximum time value of P' for which 
the slope of P'QR is less than or equal to criocai- The 
maximum Tiocai exists because the set of allowed val- 
ues of T where P' — {p, T) is closed and therefore 
compact. To compute Tiocai we substitute aiocai in the 
condition for causality of P'QR. 

Unlike the ID x Time case, it is not clear that T^^p can 
be computed by ray shooting queries. In 2D x Time, we 
need an oracle to determine which among several right 
circular cones is intersected first by a triangle P'Q'R 
when the vertex P of APQR is lifted to P' = (p, T) 
while also lifting Q to Q' = {q,t{q) -t-Tmin). However, 
just as for the ID x Time case, we can approximate 
Tiup up to any given numerical accuracy by perform- 
ing a binary search in the interval [ti{p) + Tinin, Tiocai] 
which we know contains T^^p. Therefore, the eventual 



height of the tentpole PP' is at least max{Tmin, T^up — 
■q} where ?7 > is the desired numerical accuracy. 

We thus have the following theorem. 

Theorem 7. Given a triangulation M of a bounded 
planar space domain where Wmin is the minimum 
width of a simplex of M and Omin is the mini- 
mum slope anywhere in M x [0,oo), for every e 
such that < e < ^ our algorithm constructs 
a simplicial mesh of M x [0, T] consisting of at 
most j" dzam(M) g-^m A(M) ^1 gp^jggjjj^g gZements for ev- 
ery real T > 0. 

Proof. By Theorems |4] and [5] it follows that the height 
of each tentpole constructed by the algorithm is at 
least Tmin = ecminUimin. By Theorem |6l after con- 
structing at most k < j" diamCflf) CTmax ^1 patches, the 
entire front t^ is past the target time T. Since each 
patch consists of at most A(Af) elements, the theorem 
follows. □ 

4. CONCLUSION 

We have shown how to extend the Tent Pitcher al- 
gorithm for planar and linear spatial domains to the 
case of changing wavespeeds. Our expressions for the 
causality and progress constraints that apply at each 
step make explicit the dependence on the slope of the 
cone of influence most constraining the progress at 
that step. This dependence is not explicit in the for- 
mula of Erickson et al. because they assume without 
loss of generality that the slope is 1 everywhere in 
spacetime. For the constant wavespeed case, the algo- 
rithm in this paper is an alternative to the algorithm 
due to Erickson et al. with potentially weaker progress 
constraints. We can view the algorithm of Erickson et 
al. as looking one step ahead in the sense that the 
progress constraint at step i guarantees that the front 
constructed in step i -I- 1 is causal. Our algorithm 
can be viewed as looking one step even further — our 
progress constraint at step i guarantees that the front 
constructed in step i -f 2 is causal. In a relatively 
straightforward manner, we can generalize this idea to 
looking at step i to the front in step i + h where h is 
a horizon parameter that can be chosen adaptively by 
the algorithm. It needs to be investigated whether the 
extra complexity of the algorithm for h > 2 is justified 
by a more efficient meshing algorithm overall. 

We have preliminary experimental results in ID x Time 
and a prototype with simulated physics in 2D x Time; 
more substantial empirical study is required and we 
expect to report results of such a study soon. One of 
the objectives of the study will be to explore different 
heuristics to choose which local minimum vertex to 
pitch at every step. Some heuristics, such as pitching 




Figure 5: An unstructured triangular spacetime mesh 
over a ID uniform space mesh. The space dimension 
is horizontal and time increases upwards. The slope at 
any point in spacetime is one of three distinct values: 
the minimum slope occurs in a band around the diagonal 
where the tentpoles are shortest; beyond a certain time 
value, the maximum slope occurs everywhere. 

the local minimum with the minimum slope (highest 
wavespeed), perform better than others. We have an 
extension to the current algorithm that allows pitch- 
ing at any vertex, not necessarily a local minimum. 
However, the extended algorithm is more complicated 
and it is not clear if the expected gains will be worth 
the extra computation time. 

Figures [5] and |6] illustrate spacetime meshes con- 
structed by our prototype implementation over ID and 
2D space meshes respectively. The ID x Time space- 
time mesh was constructed by pitching an indepen- 
dent set of local minima in non- increasing order of 
wavespeed. In other words, the algorithm preferred 
to pitch local minima adjacent to points on the front 
where the wavespeed was maximum (slope was mini- 
mum) . The 2D X Time mesh was constructed by pitch- 
ing a global minimum at every step. In either example, 
many more spacetime elements would be required to 
mesh the same volume if the height of every tentpole 
were constrained by the globally minimum slope. 

In higher dimensions, we have a theorem identical to 
Theorem [7] when every dihedral angle of every simplex 
is non-obtuse. We anticipate soon an analogous the- 
orem for arbitrary dimensional space domains in the 
presence of obtuse angles. 




Figure 6: An unstructured tetrahedral spacetime mesh 
over a triangulated uniform 2D grid. Time increases up- 
wards. The slope at any point in spacetime is one of two 
distinct values: the minimum slope occurs inside a circu- 
lar cone where the tentpoles are shortest, the maximum 
slope occurs everywhere else. 

Our algorithm can be modified to handle asymmet- 
ric cones, such as due to wave propagation through 
anisotropic media. In the presence of anisotropy, the 
most limiting cone constraint can be nonlocal. 

In a recent paper, Abedi et al. [I] extend TentPitcher 
to support another kind of adaptivity, where the size 
of the spacetime elements is adapted to a posteriori 
estimates of the numerical error. Abedi et al. apply 
hierarchical refinement and coarsening of the underly- 
ing one- or two-dimensional space mesh to adapt the 
spatial size of future spacetime elements. They extend 
the progress constraints of Erickson et al. to antici- 
pate future refinement and coarsening both of which 
change the shape of the elements on the front. The 
outstanding problem that we plan to consider next is 
to combine adaptivity to changing wavespeeds with 
refinement and coarsening for the case of planar space 
domains. It is quite straightforward to combine the 
progress constraints in this paper with those of Abedi 
et al. to support refinement in the presence of chang- 
ing wavespeeds. Coarsening can be done safely if each 
triangle after coarsening satisfies progress constraint 
[cTmin]- When coarsening is possible only under such 
strict constraints, we need to carefully prioritize each 
coarsening step so that the front is only as refined as 
necessary and not much more. 

Our research group is also implementing a parallel ver- 
sion of Tent Pitcher to run on multiple processors. The 
nonlocal nature of the constraints pose significant chal- 
lenges in the parallel setting. 

In many problems, the geometry of the space domain 
changes over time. There may also be internal bound- 
aries between different parts of the domain, e.g., sep- 



arating two distinct materials with different physical 
properties, and these internal boundaries may evolve 
over time. We would like to handle moving boundaries 
both internal and external. 
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